#include <iostream>
#include <fstream>
#include <iomanip>
#include "solver.h"

using namespace std;

Solver::Solver(const Config &cfg, TimeDiscr &timeAdvance, const Fluid &flows,
               const VarArray<Coord_t, 2, 0> &coordIn,
               const SpaceDiscr &spaceDiscrIn)
    : coord{coordIn}, spaceDiscrete{spaceDiscrIn},
      resReduct{cfg.GetResReduct()}, convTol{cfg.GetConvTol()},
      refMomentX{cfg.GetRefMomentX()}, refMomentY{cfg.GetRefMomentY()},
      refLength{cfg.GetRefLength()}, convgFname{cfg.GetConvgFname()},
      flowSolFname{cfg.GetFlowSolFname()}, surfFname{cfg.GetSurfFlowFname()},
      extIter{cfg.GetExtIter()},
      wrtSolFreq{cfg.GetWrtSolFreq()}, wrtConvFreq{cfg.GetWrtConvFreq()},
      timeDiscrete{timeAdvance}, flow{flows}, convFile{cfg.GetConvgFname()}

{
  if (!convFile)
  {
    cerr << "open convengence file err" << endl;
    exit(EXIT_FAILURE);
  }

  convFile << "TITLE = \" visualization of the convergence\" \n"
           << "VARIABLES = \"iter\", \"drho/drho1\" \n";
}

void Solver::FlowSolverEuler(void)
{
  int iter{0};
  real_t drho, drhoMax, drho1;
  Coord_t drhoMaxPos;

  ++iter;
  timeDiscrete.SolveOneStepEuler();
  timeDiscrete.DensityChange2Norm(drho1, drhoMax, drhoMaxPos);
  cout << drho1 << endl;
  drho1 += 1.0e-32;
  // exit(EXIT_FAILURE);

  do
  {
    ++iter;
    cout << iter << endl;
    timeDiscrete.SolveOneStepEuler();
    timeDiscrete.DensityChange2Norm(drho, drhoMax, drhoMaxPos);

    cout << log10(drho / drho1) << "\n";
    // if (iter % wrtSolFreq == 0)
    //   WriteFlowASSCII(iter);
    // cout << "okay7" << endl;
  } while ((iter < extIter) &&
           //  (convTol < log10(drho)) &&
           (-log10(drho / drho1) < resReduct));
}

void Solver::FlowSolverNS(void)
{
  int iter{0};
  real_t drho, drhoMax, drho1;
  Coord_t drhoMaxPos;

  ++iter;
  timeDiscrete.SolveOneStepNS();
  timeDiscrete.DensityChange2Norm(drho1, drhoMax, drhoMaxPos);
  cout << drho1 << endl;
  drho1 += 1.0e-32;
  // exit(EXIT_FAILURE);

  do
  {
    ++iter;
    timeDiscrete.SolveOneStepNS();
    timeDiscrete.DensityChange2Norm(drho, drhoMax, drhoMaxPos);

    if (iter % wrtConvFreq == 0)
    {
      cout << iter << "      ";
      cout << log10(drho / drho1) << endl;
    }

    convFile << iter << " " << log10(drho / drho1) << "\n";

    if (iter % wrtSolFreq == 0)
    {
      WriteFlowASSCIIPoint(iter);
      convFile.flush();
    }
    // cout << "okay7" << endl;
  } while ((iter < extIter) &&
           //  (convTol < log10(drho)) &&
           (-log10(drho / drho1) < resReduct));

  convFile.flush();
  convFile.close();
  WriteFlowASSCIIPoint(iter);
}

void Solver::WriteFlowASSCIIBlock(int iter)
{
  std::string strIter{"_" + to_string(iter) + ".dat"};
  std::ofstream flowFile{flowSolFname + strIter}, surfFile{surfFname + strIter};

  const VarArray<ConsVar_t, 2, 0> &cv{flow.GetCV()};
  const VarArray<DeptVarNS_t, 2, 0> &dv{flow.GetdvNS()};
  const VarArray<real_t, 2, 0> &p{flow.GetPressure()};
  // const VarArray<Coord_t, 2, 0> &coord{timeDiscrete.GetCoord()};

  const int iL{cv.Dim(0) - 2}, jL{cv.Dim(1) - 2};
  const int I2{iL - 1}, J2{jL - 1};
  flowFile << "TITLE = \" visualization of the 2Dflow\" \n"
           << "VARIABLES = \"x\", \"y\", \"rho\", \"u\", "
           << "\"v\", \"p\", \"T\", \"Mach\", \"lamViscous\" \n";
  flowFile << "ZONE  I=" << iL - 1
           << ", J=" << jL - 1
           << ", DATAPACKING=BLOCK, VARLOCATION=([3,4,5,6,7,8,9]=CELLCENTERED)\n";
  flowFile << std::scientific << std::setprecision(16);
  surfFile << std::scientific << std::setprecision(16);
  for (int j = 2; j <= jL; ++j)
  {
    for (int i = 2; i <= iL; ++i)
    {
      flowFile << coord(i, j).x << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;
  // 输出为 point 好些
  for (int j = 2; j <= jL; ++j)
  {
    for (int i = 2; i <= iL; ++i)
    {
      flowFile << coord(i, j).y << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << cv(i, j).rho << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << cv(i, j).rhou / cv(i, j).rho << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << cv(i, j).rhov / cv(i, j).rho << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << p(i, j) << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << dv(i, j).T << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  double u, v, mach;
  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      u = cv(i, j).rhou / cv(i, j).rho;
      v = cv(i, j).rhov / cv(i, j).rho;
      mach = sqrt(u * u + v * v) / dv(i, j).c;
      flowFile << mach << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;

  for (int j = 2; j <= J2; ++j)
  {
    for (int i = 2; i <= I2; ++i)
    {
      flowFile << dv(i, j).mu << " ";
      if (i % 10 == 0)
        flowFile << "\n";
    }
  }
  flowFile << endl;
  flowFile.close();

  surfFile << "TITLE = \" visualization of the surface\" \n"
           << "VARIABLES = \"x\", \"y\", \"rho\", \"u\", "
           << "\"v\", \"p\", \"T\", \"Mach\", \"lamViscous\", \"cf\", \"cp\" \n";
  surfFile << "ZONE  I=" << iL - 1
           << ", J=" << jL - 1
           << ", DATAPACKING=BLOCK, VARLOCATION=([3,4,5,6,7,8,9,10]=CELLCENTERED)\n";

  surfFile.close();
}

void Solver::WriteFlowASSCIIPoint(int iter)
{
  std::string strIter{"_" + to_string(iter) + ".dat"};
  std::ofstream flowFile{flowSolFname + strIter}, surfFile{surfFname + strIter};

  const VarArray<ConsVar_t, 2, 0> &cv{flow.GetCV()};
  const VarArray<DeptVarNS_t, 2, 0> &dv{flow.GetdvNS()};
  const VarArray<real_t, 2, 0> &p{flow.GetPressure()};
  const real_t qInf{flow.GetDynamicPress()}, pInf{flow.GetpInf()};

  const std::vector<BoundDescrp_t> &Bdescrp = spaceDiscrete.GetBoundaryDiscrpAll();
  const VarArray<Grad_t, 2, 0> &gradFi = spaceDiscrete.GetGradientFaceI();
  const VarArray<Grad_t, 2, 0> &gradFj = spaceDiscrete.GetGradientFaceJ();

  const VarArray<Area_t, 2, 0> &si{spaceDiscrete.GetSI()};
  const VarArray<Area_t, 2, 0> &sj{spaceDiscrete.GetSJ()};

  // const VarArray<Coord_t, 2, 0> &coord{timeDiscrete.GetCoord()};

  const int iL{cv.Dim(0) - 2}, jL{cv.Dim(1) - 2};
  const int I2{iL - 1}, J2{jL - 1};
  flowFile << "TITLE = \" visualization of the 2Dflow\" \n"
           << "VARIABLES = \"x\", \"y\", \"rho\", \"u\", "
           << "\"v\", \"p\", \"T\", \"Mach\", \"lamViscous\" \n";
  flowFile << "ZONE  I=" << I2
           << ", J=" << J2
           << ", DATAPACKING=POINT\n";
  flowFile << std::scientific << std::setprecision(16);

  real_t x, y, u, v, c, mach, rhoav, pav, tav, muav;
  for (int j = 2; j <= jL; ++j)
  {
    for (int i = 2; i <= iL; ++i)
    {
      x = 0.25 * (coord(i, j).x + coord(i - 1, j).x + coord(i, j - 1).x + coord(i - 1, j - 1).x);
      y = 0.25 * (coord(i, j).y + coord(i - 1, j).y + coord(i, j - 1).y + coord(i - 1, j - 1).y);
      u = 0.25 * (cv(i, j).rhou / cv(i, j).rho + cv(i - 1, j).rhou / cv(i - 1, j).rho +
                  cv(i, j - 1).rhou / cv(i, j - 1).rho + cv(i - 1, j - 1).rhou / cv(i - 1, j).rho);
      v = 0.25 * (cv(i, j).rhov / cv(i, j).rho + cv(i - 1, j).rhov / cv(i - 1, j).rho +
                  cv(i, j - 1).rhov / cv(i, j - 1).rho + cv(i - 1, j - 1).rhov / cv(i - 1, j).rho);
      c = 0.25 * (dv(i, j).c + dv(i - 1, j).c + dv(i, j - 1).c + dv(i - 1, j - 1).c);
      mach = sqrt(u * u + v * v) / c;
      rhoav = 0.25 * (cv(i, j).rho + cv(i - 1, j).rho + cv(i, j - 1).rho + cv(i - 1, j - 1).rho);
      pav = 0.25 * (p(i, j) + p(i - 1, j) + p(i, j - 1) + p(i - 1, j - 1));
      tav = 0.25 * (dv(i, j).T + dv(i - 1, j).T + dv(i, j - 1).T + dv(i - 1, j - 1).T);
      muav = 0.25 * (dv(i, j).mu + dv(i - 1, j).mu + dv(i, j - 1).mu + dv(i - 1, j - 1).mu);
      flowFile << x << " "
               << y << " "
               << rhoav << " "
               << u << " "
               << v << " "
               << pav << " "
               << tav << " "
               << mach << " "
               << muav << "\n";
    }
  }
  flowFile.flush();
  flowFile.close();

  surfFile << std::scientific << std::setprecision(16);
  surfFile << "TITLE = \" visualization of the surface multi zone\" \n"
           << "VARIABLES = \"x\", \"y\", \"rho\", \"u\", "
           << "\"v\", \"p\", \"T\", \"lamViscous\", \"cf\", \"cp\" \n";

  int nbound = Bdescrp.size();
  int Bncell, beg, end, count{0}, ins1, dum1, IndexGeo;
  std::string title{"wall "};
  int iins1, jins1, idum1, jdum1, igeo1, jgeo1, igeo2, jgeo2;
  real_t rho, press, T, mu, sx, sy, sxn, syn, ds;
  real_t grdnx, grdny, grdnn, dvdnx, dvdny, dvdna, cf, cp;

  for (int i = 0; i < nbound; ++i)
  {
    // 如果是壁面就输出，其他就过
    if (Bdescrp[i].BType == Bound_t::Wall)
    {
      Transij trans(Bdescrp[i].isIface);
      ++count;
      beg = Bdescrp[i].beg;
      end = Bdescrp[i].end;
      Bncell = end - beg + 1;
      ins1 = Bdescrp[i].ins1;
      dum1 = Bdescrp[i].dum1;
      // 面和点的索引与单元有差异
      IndexGeo = (Bdescrp[i].Bn < 2) ? ins1 : dum1;

      title += to_string(count);
      surfFile << "ZONE "
               << "T=\"" << title << "\", "
               << "I=" << Bncell << ", DATAPACKING=POINT\n";

      const VarArray<Area_t, 2, 0> &S{Bdescrp[i].isIface ? si : sj};
      const VarArray<Grad_t, 2, 0> &grad{Bdescrp[i].isIface ? gradFi : gradFj};
      real_t sgn{Bdescrp[i].Bn < 2 ? -1.0 : 1.0};

      for (int j = beg; j <= end; ++j)
      {
        trans(j, IndexGeo, igeo1, jgeo1);
        trans(j + 1, IndexGeo, igeo2, jgeo2);
        x = 0.5 * (coord(igeo1, jgeo1).x + coord(igeo2, jgeo2).x);
        y = 0.5 * (coord(igeo1, jgeo1).y + coord(igeo2, jgeo2).y);
        sx = sgn * S(igeo1, jgeo1).sx;
        sy = sgn * S(igeo1, jgeo1).sy;
        ds = sqrt(sx * sx + sy * sy);
        sxn = sx / ds;
        syn = sy / ds;

        // ! 速度在面法线方向的梯度
        grdnx = grad(igeo1, jgeo1).ux * sxn + grad(igeo1, jgeo1).uy * syn;
        grdny = grad(igeo1, jgeo1).vx * sxn + grad(igeo1, jgeo1).vy * syn;
        // ! 垂直于壁面分量速度的梯度
        grdnn = grdnx * sxn + grdny * syn;
        // ! 切向速度沿壁面法线方向的梯度
        dvdnx = grdnx - grdnn * sxn;
        dvdny = grdny - grdnn * syn;
        // 这就只有正的了
        dvdna = sqrt(dvdnx * dvdnx + dvdny * dvdny);

        trans(j, dum1, idum1, jdum1);
        trans(j, ins1, iins1, jins1);
        rho = 0.5 * (cv(idum1, jdum1).rho + cv(iins1, jins1).rho);
        u = 0.5 * (cv(idum1, jdum1).rhou / cv(idum1, jdum1).rho +
                   cv(iins1, jins1).rhou / cv(iins1, jins1).rho);
        v = 0.5 * (cv(idum1, jdum1).rhov / cv(idum1, jdum1).rhov +
                   cv(iins1, jins1).rhov / cv(iins1, jins1).rhov);
        press = 0.5 * (p(idum1, jdum1) + p(iins1, jins1));
        T = 0.5 * (dv(idum1, jdum1).T + dv(iins1, jins1).T);
        mu = 0.5 * (dv(idum1, jdum1).mu + dv(iins1, jins1).mu);
        // cf = tau / (0.5*rhoInf*velInf*velInf)
        // 没有定义方向也没有定义正负号  注意
        cf = 2.0 * mu * dvdna / qInf;
        cp = (press - pInf) / (qInf);

        surfFile << x << " " << y << " " << rho << " "
                 << u << " " << v << " " << press << " "
                 << T << " " << mu << " " << cf << " "
                 << cp << "\n";
      }
    }
  }
  surfFile.flush();
  surfFile.close();
}
